function [min_male, max_male, min_female, max_female, problem] = get_ricebs(couple,singlemale,singlefemale,leisurem,leisuref,hworkm,hworkf,wagem,wagef,q,Q,nonlabor,mutually_acceptable,divorce_costsM,divorce_costsF,divorce_costsMF,hcat,text_solver)

% ================
% define constants
% ================

T = 112;                                              % total potential working hours
H = length(leisurem);                                 % total number of households in the data
alpha = 0.4;                                          % bound for non-labor income

% =========================
% define decision variables
% =========================
yalmip('clear')
qM = sdpvar(H,1);                                     % private consumption of Hicksian good male
qF = sdpvar(H,1);                                     % private consumption of Hicksian good female

pM = sdpvar(H,H,'full');                              % personalized public good price male
pF = sdpvar(H,H,'full');                              % personalized public good price female

pmhworkM = sdpvar(H,H,'full');                        % personalized male's hwork price male
pmhworkF = sdpvar(H,H,'full');                        % personalized male's hwork price female

pfhworkM = sdpvar(H,H,'full');                        % personalized female's hwork price male
pfhworkF = sdpvar(H,H,'full');                        % personalized female's hwork price female

wageM = sdpvar(H,1);                                  % shadow wage male
wageF = sdpvar(H,1);                                  % shadow wage female

nl_incomeM = sdpvar(H,1);                             % non-labor income male
nl_incomeF = sdpvar(H,1);                             % non-labor income female

sharem = sdpvar(H,1);                                 % consumption based RICEB male
sharef = sdpvar(H,1);                                 % consumption based RICEB female


%%
% ==================
% define constraints
% ==================

Constraints = [];

% restriction on individual private consumption
for i=1:H    
    Constraints = [Constraints, q(i) == qM(i) + qF(i), qM(i) >=0, qF(i) >= 0];
    
    if singlefemale(i) == 1
        Constraints = [Constraints,qM(i) == 0];
    elseif singlemale(i) == 1
        Constraints = [Constraints,qF(i) == 0];
    end
    
end


% shadow price should be non-negative and equal to observed wage if the individual is employed
for i=1:H
    
    Constraints = [Constraints,wageF(i) >= 0, wageM(i) >= 0];
    
    if ~isnan(wagem(i))
        Constraints = [Constraints, wageM(i) == wagem(i)];
    end
    if ~isnan(wagef(i))
        Constraints = [Constraints, wageF(i) == wagef(i)];
    end
    
end


% restriction on public prices
for i=1:H
    for j=1:H
        
        Constraints = [Constraints, 1 == pM(i,j) + pF(i,j), pM(i,j) >=0, pF(i,j) >= 0];
        Constraints = [Constraints, wageM(i) == pmhworkM(i,j) + pmhworkF(i,j), pmhworkM(i,j) >=0, pmhworkF(i,j) >= 0];
        Constraints = [Constraints, wageF(j) == pfhworkM(i,j) + pfhworkF(i,j), pfhworkM(i,j) >=0, pfhworkF(i,j) >= 0];
        
        if (couple(i) == 1 || singlemale(i) == 1) && singlemale(j) == 1
            Constraints = [Constraints, 1 == pM(i,j), wageM(i) == pmhworkM(i,j), wageF(j) == pfhworkM(i,j)];
        end
        
        if (couple(j) == 1 || singlefemale(j) == 1) && singlefemale(i) == 1
            Constraints = [Constraints, 1 == pF(i,j), wageM(i) == pmhworkF(i,j), wageF(j) == pfhworkF(i,j) ];
        end
    end
    
end


% feasibility restrictions on the individual nonlabor income (for each matched pair)
% individual nonlabor incomes after divorce must lie between 40% and 60% of total nonlabor income under marriage
for i=1:H
    
    Constraints = [Constraints, nonlabor(i) == nl_incomeM(i) + nl_incomeF(i)];
    
    if couple(i) == 1
        if nonlabor(i)>=0
            Constraints = [Constraints, alpha*nonlabor(i) <= nl_incomeM(i), nl_incomeM(i) <= (1-alpha)*nonlabor(i), alpha*nonlabor(i) <= nl_incomeF(i), nl_incomeF(i) <= (1-alpha)*nonlabor(i)];
        else
            Constraints = [Constraints,(1-alpha)*nonlabor(i) <= nl_incomeM(i), nl_incomeM(i) <= alpha*nonlabor(i), (1-alpha)*nonlabor(i) <= nl_incomeF(i), nl_incomeF(i) <= alpha*nonlabor(i)];
        end
    elseif singlefemale(i) == 1
        Constraints = [Constraints, nonlabor(i) == nl_incomeF(i)];
    elseif singlemale(i) == 1
        Constraints = [Constraints, nonlabor(i) == nl_incomeM(i)];
    end
    
end


% individual rationality 
for i = 1:H
        if couple(i) == 1        
            Constraints = [Constraints, wageM(i)*T + nl_incomeM(i) - divorce_costsM(i) <= qM(i) + wageM(i)*leisurem(i) + Q(i) + wageM(i)*hworkm(i) + wageF(i)*hworkf(i)];
            Constraints = [Constraints, wageF(i)*T + nl_incomeF(i) - divorce_costsF(i) <= qF(i) + wageF(i)*leisuref(i) + Q(i) + wageM(i)*hworkm(i) + wageF(i)*hworkf(i)];        
        end        
end
    

% no blocking pairs
for i = 1:H
    for j = 1:H        
        if mutually_acceptable(i,j) == 1
            Constraints = [Constraints, (wageM(i) + wageF(j))*T + nl_incomeM(i) + nl_incomeF(j) - divorce_costsMF(i,j) <= ...
                qM(i) + qF(j) + wageM(i)*leisurem(i) + wageF(j)*leisuref(j) + pM(i,j)*Q(i) + pF(i,j)*Q(j) + ...
                pmhworkM(i,j)*hworkm(i) + pmhworkF(i,j)*hworkm(j) + pfhworkM(i,j)*hworkf(i) + pfhworkF(i,j)*hworkf(j)];
        end
        
    end
end


% define RICEBs
for i = 1:H
    if couple(i) == 1
        Constraints = [Constraints, sharem(i) == (qM(i) + Q(i))/(q(i) + Q(i)),...
            sharef(i) == (qF(i) + Q(i))/(q(i) + Q(i))];
    elseif singlemale(i) == 1
        Constraints = [Constraints, sharem(i) == 1];
    elseif singlefemale(i) == 1
        Constraints = [Constraints, sharef(i) == 1];
    end
end



%%
% =================
% Solve the problem
% =================

min_male = zeros(6,1);
max_male = zeros(6,1);
min_female = zeros(6,1);
max_female = zeros(6,1);
problem = 0;

options = sdpsettings('verbose',0,'solver',text_solver);

% males
j = 1;
for empym = 1:3
    for edum = 1:2

        k1 = empym*1000 + edum*100 +  1*10 + 1;
        k2 = empym*1000 + edum*100 +  1*10 + 2;        
        k3 = empym*1000 + edum*100 +  2*10 + 1;
        k4 = empym*1000 + edum*100 +  2*10 + 2;
        k5 = empym*1000 + edum*100 +  3*10 + 1;
        k6 = empym*1000 + edum*100 +  3*10 + 2;

        index = (couple == 1 & ismember(hcat,[k1,k2,k3,k4,k5,k6]));

        if sum(index) > 0 && problem == 0
            obj1 = sum(sharem.*index)/sum(index);

            % minimize males' share
            sol = optimize(Constraints,obj1,options);
            if sol.problem ~= 0
                sol.info;
                yalmiperror(sol.problem);
                problem = 1;
            else
                min_male(j) = value(obj1);
            end

            % maximize males' share
            sol = optimize(Constraints,-obj1,options);
            if sol.problem ~= 0
                sol.info;
                yalmiperror(sol.problem);
                problem = 1;
            else
                max_male(j) = value(obj1);
            end
        else
            min_male(j) = NaN;
            max_male(j) = NaN;
        end
        j = j + 1;
    end
end




% females
j = 1;
for empyf = 1:3
    for eduf = 1:2

        l1 = 1*1000 + 1*100 + empyf*10 + eduf;
        l2 = 1*1000 + 2*100 + empyf*10 + eduf;
        l3 = 2*1000 + 1*100 + empyf*10 + eduf;
        l4 = 2*1000 + 2*100 + empyf*10 + eduf;
        l5 = 3*1000 + 1*100 + empyf*10 + eduf;
        l6 = 3*1000 + 2*100 + empyf*10 + eduf;

        index = (couple == 1 & ismember(hcat,[l1,l2,l3,l4,l5,l6]));

        if sum(index) > 0 && problem == 0

            obj2 = sum(sharef.*index)/sum(index);

            % minimize females' share
            sol = optimize(Constraints,obj2,options);
            if sol.problem ~= 0
                sol.info;
                yalmiperror(sol.problem);
                problem = 1;
            else
                min_female(j) = value(obj2);
            end

            % maximize females' share
            sol = optimize(Constraints,-obj2,options);
            if sol.problem ~= 0
                sol.info;
                yalmiperror(sol.problem);
                problem = 1;
            else
                max_female(j) = value(obj2);
            end

        else
            min_female(j) = NaN;
            max_female(j) = NaN;
        end
        j = j + 1;
    end
end

return


